1. Field of the Invention
The invention relates to a method for estimating the frequency of a time signal through the use of a Discrete Fourier Transform (DFT) and interpolation between sample points of the DFT spectrum.
2. Description of the Related Art
Examples of methods for estimating frequency utilizing a Fourier Transform are used, for example, in FM CW (Frequency Modulated, Continuous Wave) Radar systems. The use of FM CW principals for level measurement is described by Dr. J. Olto in xe2x80x9cMickrowellensensor zur Fxc3xcllstandsmessungxe2x80x9d (Microwave sensor for level measurement), Sensoren-Technologie und Anwendungen, VDI Berichter 939, 1992, pages 95-100, and in xe2x80x9cMikrowellen messen Fxc3xcllstxc3xa4ndexe2x80x9d in Design and Electronik-Sensortechnik, May 1997, issue 10, pages 40-44. The same author describes digital methods for frequency determination of single and multi-frequency signals utilizing various methods of interpolation of the Fourier spectrums. By using interpolation a vastly more accurate distance measurement made through frequency estimations can result. Shown are for example interpolation using averaging or by parabolic approximation.
A process for estimating the frequency of a time signal using Fourier Transform and the interpolation between support points of the discrete Fourier transformation using a Hamming windowing is described in xe2x80x9cHighly Accurate Frequency Interpolation of Apodized FFT Magnitude-Mode Spectra,xe2x80x9d by Goto in Applied Spectroscopy, Vol. 52, Nr.1, 1998, page 134 et seq.
Evaluation methods for precision distance measurement with FM CW systems and their use in the microwave field are described by Stolle, Heuerman and Schiek in tm-Technishes Messen 62 (2/95), pages 66-73.
Methods for accurate frequency estimation of a Fourier transformed time signal through interpolation are described by Jain, Collins and Davis in xe2x80x9cHigh-Accuracy Analog Measurements via Interpolated FFTxe2x80x9d IEEE Vol. IM-28, No. 2, June 1979, pages 1213-122 and by Grandke in xe2x80x9cInterpolation Algorithms for Discrete Fourier Transforms of Weighted Signals,xe2x80x9d IEEE Vol. IM-32, No. 2, June 1983, pages 350-355.
Further methods of frequency estimation can be found in xe2x80x9cCalculation of Narrow-Band Spectra by Direct Decimationxe2x80x9d by Liu and Mintzer in IEEE Transactions, Vol. ASSP-26, No. 6, December 1978, pages 529 through 534; and in xe2x80x9cSome Aspects of the Zoom Transform,xe2x80x9d by Yip in IEEE Transactions on Computers, Vol. C-25 No. 3, March 1976, pages 297 through 296.
The precise frequency measurement is carried out by the known Fourier-spectrum synthesis methods primarily accomplished through interpolation formula for the signal filtering through rectangular- or Hanning-windowing functions that are describable as simple equations. The Hanning interpolation formula can also be described as a signal with Hamming windowing under the assumption of an infinite expansion of a complex time signal and accepting the systematic errors of that formula.
The invention is based upon providing a method of frequency estimation that is advantageous with regards to frequency accuracy and processing expense.
The method of the invention takes advantage of the well known advantageous spectral characteristics of the Hamming window function and in particular allows an increase of the precision of a frequency estimate using a single detection with less complexity by not requiring insertion of the numerical method for the determination of the interpolation size. In particular, the invention may be represented as a closed end solution of a systematically correct equation for infinitely long real mono-tone signals as an approximation of real signals of finite length.
With the method of the invention, there are no further required approximations for the determination of the interpolation size, for instance through iterative, numerical evaluation of an interpolation rule. The off-line construction of interpolation weights, whose tabular arrangement is certainly possible and advantageous, the closed form solution algorithm will preferably be carried out as an online calculation with a default. By default, a tabular evaluation is more advantageously the maximum deviation of the stored values dependent upon the solution determined from the ratios a from adjacent maximum to high-maximum and is preferably smaller than the m-th part of the maximum values of the interpolation size with m being the number of the solution increment of the defined range of xcex1.
The invention is further illustrated below by derivation and presentation of a preferred algorithm.
For the sake of simplicity, a single frequency timing signal of frequency f0 will be examined as a discrete time signal, which is inside of a limited time window of length Nxc2x7TA in the form of N discrete sample values s(k) with a distance of the sampling period TA.
s(k)=Vxc2x7sin(2xcfx80ƒ0xc2x7kxc2x7TA+"PHgr")0xe2x89xa6kxe2x89xa6Nxe2x88x921
with V defined to be the amplitude and "PHgr" defined to by the initial phase (k=0). The generalized cosine window       w    ⁡          (      i      )        =      a    -                            (                      1            -            a                    )                ·        cos            ⁢              xe2x80x83            ⁢              (                  2          ⁢                      π            ·                          i              N                                      )            
xe2x80x83of the Hamming window function with a=0.54 has the Fourier transform function:       W    ⁡          (      k      )        =      N    ·          (                                                  -                                                1                  -                  a                                2                                      ·            δ                    ⁢                      xe2x80x83                    ⁢                      (                          k              +              1                        )                          -                                                            1                =                a                            2                        ·            δ                    ⁢                      xe2x80x83                    ⁢                      (                          k              -              1                        )                              )      
xe2x80x83with k as the discrete sample (bin) of the spectrum and xcex4 being the dirac-delta function.
Through the convolution of the Fourier-transform of the Hamming window function W(k) with the timing signal s(k) and application of the discrete Fourier Transform, one obtains the DFT-spectrum Sw(i) which is the Hamming window weighted time signal with i as the number of the discrete sample (or bin) of the spectrum.
As a rule, the frequency f0 of the timing signal does not exactly match the frequency of the discrete frequency bin of the DFT spectrum fi=il(N*TA). For the frequency f0 which is not known a priori, many lines will result in the DFT spectrum, with a highest maximum at i=l as the highest value line and one of the adjacent (largest) neighboring maximum appears at i=l+1 or i=lxe2x88x921. The sought after frequency f0 lies in the frequency interval between sample (or bin) i=lxe2x88x921 and sample (or bin)i=l+1 at a correction distance d from the main line i=l according to f0=xcex* (N*TA) with xcex=l+d and xe2x88x921xe2x89xa6dxe2x89xa6+1, where d will also be denoted as a sub-bin.
Preferably the more exact frequency estimate based upon the correction distance calculation will be based upon the size of the highest maximum Sw(i) and the size of the adjacent maximum Sw(ixc2x11). Preferably, the ratio of these two sizes   α  =            |                        S          w                ⁡                  (                      l            ±            1                    )                    |              |                        S          w                ⁡                  (          l          )                    |      
xe2x80x83will be formed, which is in the range of 0xe2x89xa6xcex1xe2x89xa61. For the application of these ratios as an aid, the relation is more advantageously derived for the solution of the correction distance d,   α  =            |                        (                                    -              2                        +                          4              ⁢              d                        -                          2              ⁢                              d                2                                      +                          a              ·                              (                                                      -                    2                                    -                                      8                    ⁢                    d                                    +                                      4                    ⁢                                          d                      2                                                                      )                                              )                ·                  (                      d            +            1                    )                    |              |                                    -            2                    ⁢          d                +                  a          ·                      (                                          4                ⁢                                  d                  2                                            -              2                        )                    ·                      (                          d              -              2                        )                              |      
xe2x80x83which can be simplified by merely making small angle approximations. According to the definition of the Hamming window, a=0.54 and the prominent relation for xcex1 transforms itself into the third degree equation
d3(0.16xcex1+0.16)+d2(xe2x88x920.32xcex1xe2x88x920.16)+d(xe2x88x921.08xcex1xe2x88x921.24)+2.16xcex1xe2x88x920.92=0
This equation advantageously represents a closed end solution for the defined region of xcex1 (0xe2x89xa6xcex1xe2x89xa61), according to which the calculation of the correction distance will be preferably realized. The introduction of an aid to simplify notation results in the closed end solution for   d  =            x      3        -                  a        2            2      
where       x    =                  -        2            ·                                    |                          p              3                        |                    ⁣                                    ·              cos                        ⁢                          xe2x80x83                        ⁢                          (                                                ϕ                  3                                +                                  π                  3                                            )                                                ϕ    =          arccos      ⁢              xe2x80x83            ⁢              (                              -                          q              2                                ·                                    (                                                |                  p                  ⁢                                      |                    3                                                  3                            )                                      )                  p    =                  a        1            -                        1          3                ⁢                  a          2          2                          q    =                            2          27                ⁢                  a          2          3                    -                        1          3                ⁢                  a          1                ⁢                  a          2                    +              a        0                        a      0        =                  (                              2.16            ⁢                          xe2x80x83                        ⁢            α                    -          0.92                )            c                  a      1        =                  (                                            -              1.08                        ⁢                          xe2x80x83                        ⁢            α                    -          1.24                )            c                  a      2        =                  (                                            -              0.32                        ⁢                          xe2x80x83                        ⁢            α                    -          0.16                )            c            c    =                  0.16        ⁢                  xe2x80x83                ⁢        α            +              0.16        .            
The prominent equation for d forms an analytical and definite solution for the problem of the interpolation between sample points (or bins) of the DFT spectrum with a Hamming Window function filtered time signal on the basis of the ratios of the highest maximum and its adjacent maximum.
The sequence of Sw(l) for the highest maximum and Sw(l+1) for the adjacent maximum forms the basis for the derivation of the solution for d. When the location of the adjacent maximum is i=lxe2x88x921, the accurate estimation of the frequency is xcex=l+d, merely the sign of the ascertained value for the correction distance d is inverted.